#include <stdio.h>
#include <string.h>
#include <math.h>

#include "common.h"

#define CPUTILESIZE 32

/*
 * Safe version using GCC's typeof to avoid double evaluation problem, e.g min(++a,++b)
 */
#define max(a,b) \
   ({ __typeof__ (a) _a = (a); \
       __typeof__ (b) _b = (b); \
     _a > _b ? _a : _b; })

#define min(a,b) \
   ({ __typeof__ (a) _a = (a); \
       __typeof__ (b) _b = (b); \
     _a < _b ? _a : _b; })


/*
 *  CPU version
 *  No optimization
 */
void corrcoefmatrix0(float *out, float *in, int n, int m) {

  int i,j,k;
  float a_x,b_x;
  float sum_a, sum_b, sum_a2, sum_b2, sum_ab, corrcoef;

  for (i=0; i < n; i++) {
     out[i*n+i] = 1.0;
     for (j=i+1; j < n; j++) {
        sum_a = sum_b = sum_a2 = sum_b2 = sum_ab = 0;
        for (k=0; k < m; k++) {
           a_x = in[i*m+k];
           b_x = in[j*m+k]; 
           sum_a += a_x;
           sum_b += b_x;
           sum_a2 += a_x*a_x;
           sum_b2 += b_x*b_x;
           sum_ab += a_x*b_x;
        }
        corrcoef = (m*sum_ab - sum_a*sum_b)/sqrtf((m*sum_a2-sum_a*sum_a)*(m*sum_b2-sum_b*sum_b));
        out[i*n+j] = corrcoef;
        out[j*n+i] = corrcoef;
     }
  }

}

/*
 *  CPU version
 *  Moving sum_a, sum_a2 to the outer loop
 */
void corrcoefmatrix(float *out, float *in, int n, int m) {

  int i,j,k;
  float a_x,b_x;
  float sum_a, sum_b, sum_a2, sum_b2, sum_ab, corrcoef;

  for (i=0; i < n; i++) {
     out[i*n+i] = 1.0;
     sum_a = sum_a2 = 0;
     for (k=0; k < m; k++) {
        a_x = in[i*m+k];
        sum_a += a_x;
        sum_a2 += a_x*a_x;
     }
     for (j=i+1; j < n; j++) {
        sum_b = sum_b2 = sum_ab = 0;
        for (k=0; k < m; k++) {
           a_x = in[i*m+k];
           b_x = in[j*m+k]; 
           sum_b += b_x;
           sum_b2 += b_x*b_x;
           sum_ab += a_x*b_x;
        }
        corrcoef = (m*sum_ab - sum_a*sum_b)/sqrtf((m*sum_a2-sum_a*sum_a)*(m*sum_b2-sum_b*sum_b));
        out[i*n+j] = corrcoef;
        out[j*n+i] = corrcoef;
     }
  }

}

/*
 *  CPU version with tiling to test cache
 *  
 */
void corrcoefmatrix_tiled(float *out, float *in, int n, int m) {

  int i,j,k;
  float a_x,b_x;
  float sum_a, sum_b, sum_a2, sum_b2, sum_ab, corrcoef;
  int ii,jj;
  int last_i,last_j; 
  int first_j;
  int numtiles = (n+CPUTILESIZE-1)/CPUTILESIZE;
  int numelems = 0;

  for (ii = 0; ii < n; ii += CPUTILESIZE) {
     for (jj =0; jj < n; jj += CPUTILESIZE) {
         last_i = min(ii+CPUTILESIZE,n);
         last_j = min(jj+CPUTILESIZE,n);
//numelems = 0;
//fprintf(stderr,"ii=%d-%d, jj=%d-%d\n",ii,last_i,jj,last_j);
         for (i = ii; i < last_i; i++) {
            out[i*n+i] = 1.0;      // This is repeated many times, one for each tile.
            sum_a = sum_a2 = 0;
            for (k=0; k < m; k++) {
               a_x = in[i*m+k];
               sum_a += a_x;         // This is also repeated.
               sum_a2 += a_x*a_x;    // This is also repeated.
            }
            first_j = max(jj,i+1);
            for (j=first_j; j < last_j; j++) {
               sum_b = sum_b2 = sum_ab = 0;
               for (k=0; k < m; k++) {
                  a_x = in[i*m+k];
                  b_x = in[j*m+k]; 
                  sum_b += b_x;
                  sum_b2 += b_x*b_x;
                  sum_ab += a_x*b_x;
               }
               corrcoef = (m*sum_ab - sum_a*sum_b)/sqrtf((m*sum_a2-sum_a*sum_a)*(m*sum_b2-sum_b*sum_b));
               out[i*n+j] = corrcoef;
               out[j*n+i] = corrcoef;
//numelems++;
            }
         }

//fprintf(stderr,"%d\n",numelems);
     }
  }

}
